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Abstract 

The algebraic formulation of Large N matrix mechanics recently 
developed by Halpern and Schwartz leads to a practical method of nu- 
merical computation for both action and Hamiltonian problems. The 
new technique posits a boundary condition on the planar connected 
parts X w , namely that they should decrease rapidly with increasing or- 
der. This leads to algebraic/variational schemes of computation which 
show remarkably rapid convergence in numerical tests on some many- 
matrix models. The method allows the calculation of all moments of 
the ground state, in a sequence of approximations, and excited states 
can be determined as well. There are two unexpected findings: a large 
d expansion and a new selection rule for certain types of interaction. 
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1 Introduction 



Large N matrix mechanics differs from ordinary quantum mechanics in 
that the canonical commutator 

z[p, <?] = /, (1.1) 
in the one-matrix case, is replaced by the relation 

z'M]=|0)<0| (1.2) 

where | 0) is the ground state in the reduced Hilbert space. The original 
matrix- valued coordinates <p rs , r,s = 1 . . . N are represented by the single 
operator in this reduced Hilbert space. 

The solution of the one-matrix Large N Hamiltonian problem with an 
arbitrary potential V(4>) was given some years ago H]; and only a couple of 
two-matrix problems in the action formalism have previously been solved. j|, 

§ 

The many-matrix problem involves several noncommuting operators </> m 
and their conjugate momenta. Following Halpern and Schwartz 0, this 
system is described at equal times by a symmetric free algebra which involves 
a pair (tilde and untilde) for each hermitian operator 

[4>m, 4>n] = [vTm, 7T n ] = 0, TO, 7i = 1 . . . d (1.3a) 

i[^rn,(pn] = i{-Xm,<f>n} = <W, | 0) (0 | (1.3b) 
0m | 0) = (f) m | 0), 7T m | 0) = 7T TO | 0) (1.3c) 

and the ground state energy is given by 

1 d 

Eo = N 2 (0 | - £ n m 7r m + V{4>) | 0) (1.4) 

* m=l 

where (0) refers to the set of operators {0 m }. We shall use the summation 
convention in what follows. 

In ordinary quantum mechanics systems of several interacting bodies are 
most commonly attacked from the Schrodinger equation in coordinate space, 
using the direct product basis | qi, q%, . . . qj). That approach is not available 
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in the Large N reduced Hilbert space because of the noncommutativity of 
the operators m . A basis of states in this reduced space may be written as 

\w}=(j) w \0} (1.5) 

where we use the "word" notation for ordered products of operators 

(fi w = 4> mi (fim2 ■ ■ ■ (fim n , w = m 1 m 2 ...m n , rrii = l...d (1.6) 

and we write [w] = n for the length of the word w. See App. A for a collection 
of relevant definitions and formulas. 

The new approximation technique presented in this paper lies close to the 
Heisenberg (matrix) formulation rather than the Schrodinger (wavefunction) 
formulation and makes use of the set of polynomials T w ((f>) introduced in 
Ref. 0. 

(1-P m <l> m + X{p))- X = ]T/rT w (0) (1.7a) 

w 

X(P) = X o = 0, (0 I T w (<f)) | 0) = 5 W , (1.7b) 

where the f3 m are a dummy set of (noncommuting) parameters and the num- 
bers X w were identified as the planar connected parts defined in earlier dia- 
grammatic studies. H Various properties of these X w are given in App. A, 
including their relation to the ordinary moments Z w = (0 | <p w | 0) of the 
ground state. 

The core idea of the present work is to truncate the set of these X's 

set X w = for all [w] > n (1.8) 

and solve the (now finite) set of algebraic equations, calling this the "n- 
th order approximation". Then increase n, step by step, and see whether 
the numerical results appear to converge. This is an intuitive/experimental 
approach for now, since we have no mathematical proof that this method 
should work. 

With even a small number of the X's determined, one can approximate all 
the moments of the ground state and the accuracy of these results increases 
systematically as one proceeds to higher orders of approximation. The ex- 
cited states of a Hamiltonian system are also amenable to this method. 
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The recent algebraic developments by Halpern and Schwartz || [7]] pro- 
vide a wealth of formal definitions and relations for many-matrix problems, 
unifying the study of both action and Hamiltonian systems. These start 
with the definitions of generalized creation and annihilation operators in the 
reduced Hilbert space, 

n m | 0) = iF m {<P) | 0), (0 | vr m = -i(0 | F m {<j>) (1.9a) 

B m = F m (0) + l7T m , B m | 0) = (0 | B\ n = (1.9b) 

B m B\ = E mn {4>) (1.9c) 
E mn {4>) I 0) = 2i[TX n ,F m {(j>)) | 0) (1.9d) 

which is the Interacting Cuntz Algebra. (In the case of non-interacting har- 
monic oscillators, we have E mn oc 5 mn and Eqs. (|1.9b| , |1 -9cj) reduce to the 
original Cuntz algebra.) 

In the practical work of this paper there is a basic distinction between 
the two types of problems. For action problems we start out knowing the 
functions F m (<f)) explicitly and this lets us work directly with the algebraic 
equations for the connected parts X w derived in Ref. |7| (see Sec. 5). For 
Hamiltonian problems we do not know F m (<f>) beforehand and so part of 
the method presented here involves a constructive representation of these 
operators, for which task we use the polynomials T w ((p) (see Sec. 6). 

In Sec. 2 we test the idea on a simple example: a one-matrix action 
problem and in Sec. 3 we try to give some understanding of why this method 
apparently works well. Counting of the variables in many-matrix problems 
and making use of symmetry to keep things manageable is discussed in Sec 4, 
followed in Sec. 5 by some algebraic results for a model action problem with 
d interacting matrices. The plan of attack for many-matrix Hamiltonian 
problems is set out in Sec. 6 and numerical results for a set of model potentials 
are presented in Sec. 7. We note not only the extremely rapid convergence 
found in these examples but also an unexpected selection rule. Section 8 
presents more details of this computational program; and a related method 
for calculating excited states is given in Sec. 9. Several appendices discuss 
further details and possible extensions of this work. 



3 



2 First Test: One-Matrix Action Problem 



We start with a simple problem: a one-matrix action at Large N. As given 
in Ref. f7| for the quartic action (F = 3 ), we have the following equation 
for the connected parts X n : 

x(x + i) 2 -p 2 x 2 -^p 4 = o, x = 1 £p n x n (2.1) 

1 n>0 

which leads to the recursion formula 

n~2 n-p—2 

X n = -5 n ,4 ~ X p ( 2J X q X n - p - q + 2X n _ p ), n = 4, 6, . . . . (2.2) 

2 p= 2 g=2 

For one-matrix problems we replace the word label w by n = [w]. We can 
compare this with the Schwinger-Dyson equations for the ordinary moments 
Z n = (0 I 4> n I 0), which may be written as 

n 

2^n+4 — ^ Z m Z n - m , Zq = 1 (2-3) 
m=0 

and only even n enter because of the parity symmetry in this problem. If we 
have the value of X2 = Z2 (which we know from other analysis to be (2/3) 3 / 2 ), 
then we can compute all the higher ones from these recursion formulas. Table 
1 shows some numerical results and we see that the ratio X n /Z n decreases 
fairly rapidly as n increases. 



Table 1. X n and Z n for F = <p 3 action problem 



n 


2 


4 


6 


10 


20 


X n 


.544331 


-.0925926 


.0403208 


.0143736 


-.00311591 


Z n 


.544331 


.500000 


.544331 


.816497 


3.95996 


X n j Z n 


1.00000 


-.185185 


.074074 


.017604 


.000787 



Now we want to turn this process around and calculate the value of Xi 
from the recursion formula (|2.2|) using the idea that X n should decrease 
rapidly at large n - a sort of boundary condition. That is, we consider X2 as 
an unknown parameter and then search for that value that will allow us to 
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truncate the equations ( |2.2| ) with X n+2 = 0; and then we step up the value 
of n and repeat the process. Table 2 contains the results of this computation 
and we see that the residual error at each level of approximation decreases 
quite rapidly as we increase n. 



Table 2. Compute X 2 by truncation: X n+2 = 



n+2 


4 


6 


8 


10 


20 


Approx. X 2 


.500000 


.534522 


.541429 


.543344 


.544321 


Error 


-.044331 


-.009809 


-.002902 


-.000987 


-.000010 



We view this as a sort of eigenvalue problem for the connected parts 
X n and recognize a certain similarity here with the familiar procedure for 
numerical integration of the one-dimensional Schrodinger equation in some 
given potential. While that other problem involves a continuous variable 
ijj(x) obeying a linear (differential) equation our current problem involves a 
discrete set X n obeying a nonlinear (algebraic) equation. 

3 Why Should this Method Work? 

To understand what is going on here it may help to consider the ordinary 
moments 

Z n =(O|0 n |O) = Jdqp(q)q n (3.1) 

for a one-matrix problem. These Z n , for a typical ground state density p(q), 
are a rather monotonous sequence of numbers. The infinite set of coupled 
equations for these moments (Schwinger-Dyson equations in one language) 
contains all the information about the ground state; but one would not try 
to truncate this infinite system of equations by setting the Z n equal to zero 
after some cutoff n = n*. 

(In earlier work || on moment equations for the one- and two-body 
Schrodinger equation, the asymptotic behavior of these moments as n — > oo 
was inferred from the differential equation for the wavefunction and this al- 
lowed a backward iteration procedure.) 

Now, by contrast, observe the definition of the planar connected parts, 
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again for the one-matrix problem: 

X n+l = (0 | 0T n (0) | 0) = Jdq p{q) q T n {q) (3.2) 
where the polynomials T n have the property 

(0 | T n (<j>) | 0) = n > 0. (3.3) 



Clearly, the X n are just an algebraic combination of the Z n . But Eq. (J3T 
tells us that the polynomials T n are oscillatory within the domain of inte- 
gration; and this suggests that the X n , given by (|3.2j), can be thought of 
as something like the Fourier coefficients of the density p(q). Therefore, if 
the ground state is reasonably smooth and the polynomials T n are reason- 
ably "appropriate" , then we would expect that the higher Fourier coefficients 
(the X n ) could decrease rapidly. This is the motivation to try a truncation 
scheme on the X's. 

A further advantage of the X's is that they are directly sensitive to the 
interactions in many-matrix problems. In Ref. ]7j it was shown that in many- 
matrix problems without interactions, the X w vanish if there is any mixing 
of letters in the word w. 

Once one has determined, approximately, even a small number of the X's, 
this allows one to give approximate values for all of the Z's in any one- or 
many-matrix problem by use of the general algebraic relation ( A.2j) between 
the generating functions for these two sets of parameters. 

With these encouraging results, we go on to study the problems of many 
matrices in Large N action and Hamiltonian systems. 



4 Many Matrices - Counting the Variables 

With d matrices, the number of words of length n is d n and this number 
grows very rapidly. If we have some symmetries in the action or the Hamil- 
tonian, then we can reduce the number of independent variables X w that 
we have to handle at each level of approximation. In this paper we consider 
model problems with the following invariance properties of the ground state 
|0>. 

Parity Symmetry: Change the sign of m (and 7t m ) for any m. 
Permutation Symmetry: Make any permutation among the d labels m,n, . . . . 
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In addition, there is the general invariance of the X w (as of the Trace opera- 
tion in the unreduced space) under a cyclic permutation of the letters in the 
word w. 

With these conditions, the number of independent X w 's is greatly reduced 
- to what we shall call a set of "basic words" at each level n - as shown in 
Table 3. 



Table 3. Count of d n — > basic words 



n 


d= 


=2 


d=3 


d=5 


d=9 


2 


4- 


►1 


9->l 


25-^1 


81->1 


4 


16- 


^3 


81^3 


625^3 


6561^3 


6 


64- 


+4 


729^9 


15625^9 


531441^9 


8 


256- 


+12 


6561^41 


390625^59 




10 


1024- 


^28 


59049^257 






12 


4096- 


->94 









At each level of approximation (signified by the maximum word length 
n) we shall deal with a number of basic words (the dimension D of our 
parameter space). From Table 3 we read off these dimensions: for d=2, 
D=l,4,8,20,48,. . . ; for d=3, D=l,4,13,54,. . . ; for d=5, D=l,4,13,72,. . . ; for 
d=9, D=l,4,13,. . . . The first task of the computer program is to make a 
table of all d n words at each n, identify each word with an equivalence class 
according to the symmetries described above and assign one member of each 
class as a basic word Wi, i = 1 . . . D. 

5 Many- Matrix Action Problems 

5.1 General algebraic machinery 

For action problems, we have the dual basis system of equations derived by 
Halpern and Schwartz J7|: 

Bl = G m {4>) - E mn {<j))B n , m = B m (l + X(flt)) (5.i a ) 

X(flt) = Y. X * B]W = £^ st ™ (5-lb) 
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Here, the operators B m , B^ obey the simple Cuntz algebra 

BmB^ = 5 mn (5-2) 

and the role of these operators is to generate an infinite set of coupled al- 
gebraic equations for the connected parts X w , as will be shown by example 
below. The functions G m = 2F m and E mn , defined earlier in (|1.9a| , |1.9d|) , 
are immediately known once we specify the action S. Then we shall proceed 
with the sequence of truncation approximations, generalizing the one-matrix 
example of Sec. 2. 



5.2 A model problem 

We take for our model problem here the d-matrix action 



J2 Tr([0 m ,0 n ]) 2 (5.3) 



AN 

^ JV m ,n=l 

in the unreduced Hilbert space. This gives us the reduced operators, 

G m{4>) = (0m0n0n + 4>n4>n4>m ~ ^n^m^n) (5.4a) 

E mm {<t>) = {Mn + X nn ) (5.4b) 

E m ^ n {4>) = 4>m4>n ~ 24>n4>m (5.4c) 

where we note that this S has the symmetries mentioned in the previous 
section and this leads to the simplifications X m = 0, X mn = 5 mn Xu- 
The equations (|5.1a|) now look like this 

Bl= E { 

n^m=l 

(B n B n B m X + B m (B n B n X — X n ) — 2B n B m B n X) + 
(B n XB n B m X + B m XB n B n X — 2B n XB m B n X) + 
(B n B n XB m X + B m B n XB n X — 2B n B m XB n X) + 
[B n XB n XB m X + B m XB n XB n X — 2B n X B m X B n X)} . (5.5) 
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This system of equations is equivalent to the Schwinger-Dyson set of 
equations but it is packaged to emphasize the role of the X's and it leads 
directly to our sequence of approximations. The first line of terms in ( |5.5| ) 
has only one X and its first few terms are 

(X mnn p-\- X nnnl p 2,X nmn pjBp-\- (X mnn pq r ~\~ X nnm pq r 2,X nmn pg r ) B r By B^ (5.6) 

where the usual constraint on the sum (n ^ m) is understood. The second 
and third lines have two X's and their first few terms are 

^XnXfYipBp -\- (2,XnX m pq r -\- X n pX mrl q r -\- X m pX nn q r 

^X n pX nnl q r -\- X nn pqX mr -\- X nnl pqX nr 2,X mn pqX nr ^ B r B ^B^ (5.7) 

and the fourth line, with three X's, starts off as 

(X n pX n qX mr -\- X m pX n qX nr 2,X n pX m qX nr ^ B^,BqBp. (5.8) 

Collecting the linear terms in B* gives us the equation 

1 = 2(d - 1)(X 1122 - X 1212 + X 2 n ) (5.9) 

where we have used the symmetry properties to list the basic words: (11) at 
n = 2 and (1111), (1122), (1212) at n = 4. This equation is exact and leads 
to our lowest (second) order approximation: we set all X's with word length 
greater than 2 equal to zero and we get 



X n ~l/y/2(d-l). (5.10) 

Next, we collect the cubic terms in BK For our fourth order approxima- 
tion we drop all X w 's with [w] > 4: 

X\\^(2d 2 -\- € m p -\- € mr )X m pq r 2,€ m pX m g r p '2€ mr X mr pq -\- 

(SmpSqr + S mr 5pq)[(d — l)Xn 2 2 + Xim — X mrn q q ] } + 
"^11 l^pq^mr^mp ^qr^mp^mq 2 5p r S rn q Cfjip ] 

where e pq = 1 — 5 P q. These equations are now evaluated for varying choices 
of the labels m,p, q, r, which must be paired. We find 

X mi + X1122 = form = p = q = r (5.12a) 
Xim + 3(d - 1)X 1122 - 2X 1212 + X^ = form = p^q = r (5.12b) 
2X 1122 - dX 1212 + X\ x = form = q^p = r. (5.12c) 
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The solution of this set of equations (for d 7^ 2) is 



-i^im — — -X1122 — 2^1212 — 3^ _|_ 2^ n (5.13) 

and, putting these results back into (|5.9|) , we find the fourth order approxi- 
mation for Xu 

X u ~[2(d-l)(l-^i-^)]-i. (5.14) 

For d = 2 the equations (|5.12j ) are indeterminate; but for this case a 
scaling argument leads to the conclusion that the system is not bounded. 

It was very pleasing to find, in the fourth order calculation above, that the 
number of independent equations was just equal to the number of unknowns 
and we found a unique solution. Will this circumstance continue at higher 
orders of approximation? I have no general answer. 

One should program a computer to carry the above sequence of approxi- 
mations to higher order; only algebraic work is required at each step. I have 
not done this yet, giving priority to the more difficult Hamiltonian problems, 
reported in Sec. 6. 



5.3 A large d expansion 

From the result above one is led to speculate that this truncation sequence 
of approximations may be related to a "large d" expansion. The algebraic 
calculations described above have been carried out to the sixth order, with 
9 equations in 9 unknowns, and solved in the approximation that d » 1. 
This leads to the following result: 

(X,)- 2 = 2( rf - - ^ - i§ + 0( rf -)]. (5.15) 

We do not have a systematic theory of such a large d approximation but 
the following crude attempt may be instructive. Look back at the formula 
for G m , Eq. (|5.4a| ), and replace the operator pair <p n <p n by its ground state 



average, which is X\\. This butchered G m is then 

G m ~ 2u4> m , lo = (d - l)Xn (5.16) 

which is the formula for a system of noninteracting harmonic oscillators. The 
oscillator result X mn = 6 mn /(2u) then gives immediately the leading term 



10 



in (|5.15|) . The higher order terms in 1/d are then expected to come from a 
perturbation theory expansion about this oscillator approximation. Also, if 
one looks at the computer results for the Hamiltonian problems (Sec. 7) one 
may discern a suggestion of more rapid convergence for larger values of d. 

6 Many- Matrix Hamiltonian Problems 

6.1 Choosing the model problems 

We shall study the Hamiltonians for d bosonic matrices, given in the unre- 
duced Hilbert space as 

if = ~ £ Tr{ Vm ) + NTr(V(^=)) (6.1) 

m=l V-'' 

with the following choices of the potential: 





i d 

7 E 0m 
^ m=l 


(6.2a) 


vm = 


z(E0 m ) 2 

4 m=l 


(6.2b) 


VM = 


1 d 

7 E 0m0n 
m<n=l 


(6.2c) 


V A {<j>) = 


1 d 

m<n=l 


(6.2d) 



or, if desired, any linear combination of them. The first potential, which is 
just the no n- interacting case, is used for verification of the computational 
procedure. The third and fourth potentials have "flat directions", which 
make them particularly interesting. (Will the calculations converge nicely, 
indicating a bound state, or will they not?) All four potentials have the 
symmetries (parity and permutation) described in Sec. 4. The additional 
SO(d) symmetry of V% and V4 is not used at the outset but will be noted in 
the results. 

The following subsections outline the method and further details are given 
in Sec. 8 and in Apps. A and B. 
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6.2 Construction of F m (<p) 

A central construct of our previous work j6|, [7j is the reduced operator F m (4>), 
defined in ( |1.9a|) . We will represent this quantity by a finite linear expansion 
in the polynomials T w (<p) 

F m (0) = £j& m) T«(0) (6.3) 

w 

at each level of approximation and then see how to determine the coefficients 
R. (See Subsec 8.1 for more details.). 

For any reduced operator A which depends on the 0's one has the identity 

2(0 | A{(j>)F m {(j>) | 0) = (0 | z[7f m , A{<f>)\ | 0) (6.4) 

which is proved using the definition (|1.9a|) and ( |1.3c|) . Choosing A = T w i and 
using the formulas (|A.7|) and ( |1.7b| ) this gives 

(0 | T w ,{4>)F m {4>) I 0) = ~<W,m (6.5) 

for any word w' . We impose these relations on the approximate expansion 
73[) and obtain 

£#<^ m) = |<W,m (6.6) 

where 

^={0 1^(0)^(0)10}. (6.7) 

This matrix K is numerically evaluated in terms of the X's, as detailed in 
App. B, and then we determine the expansion coefficients R from a straight- 
forward matrix inversion calculation. Of course, we make this a square (and 
positive) matrix, as detailed in Eqs. (|8.5| , |8lf ). This completes the first part 
of the fitting problem, which we would term the kinematic part since it as- 
sures that we are doing our best, at any given level of approximation, to 
represent the basic commutator algebra ( |1.3b ). 

Now we turn to the second part, which involves the dynamics of any 
particular Hamiltonian. 



12 



6.3 Minimizing the Energy 

The kinetic energy of the ground state can be expressed as 

K.E./N 2 = ^(0 | vr m vr m | 0) = ^(0 | F m F m | 0) 

= \(0\z[n m ,F m }\0) = \R^ = ±R? (6.8) 

using the methods and results of the previous subsection. 

The potential energy of the ground state is expressed directly in terms of 
the X's using flOTcp : 

(O|0fjO) = X llu + 2X 2 n (6.9a) 
(0|On|0> = ^ii23 + ^ii, rn^n (6.9b) 

(0 | <pm<Pn<Pm<Pn | 0) = X1212, m^U (6.9c) 

where we have used the specified symmetries to write these formulas in terms 
of the four basic words at the second and fourth orders. 

We program the computer to evaluate the ground state energy E = 
Eq/N 2 at the n-th order approximation with any assigned numerical val- 
ues for the quantities X w for [w] < n. The final step of this scheme is to vary 
this set of X's so as to minimize E. This procedure is without mathematical 
justification; it just seems like the natural thing to do. 

What is more, this part of the method is far from straightforward as a 
computational task because the energy E is a very nonlinear function of the 
many variables X. In Subsec. 8.2 we describe the techniques used to search 
for this minimum. The numerical results are presented next. 

7 Numerical Results 

The Tables that follow give the outputs of the computations and are designed 
to show at a glance the convergence of the approximation scheme described 
above. 

Table 4 shows the energy (E/d) calculated for the potential V 2 , for several 
values of d and at several levels of approximation and Table 5 gives the 
corresponding values of Xn = (0 | <p\ | 0). 
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Table 4. Calculated values of E/d for potential V 2 



n 


D 


d=2 


d=3 


d=5 


d=9 


2 


1 


.429 


.472 


.5408 


.6412 


4 


4 


.42672 


.47035 


.53921 


.64007 


6 


8,13 


.426672 


.4703152 


.539189 


.640058 


8 


20,54 


.42667093 


.47031461 






10 


48 


.426670885 









Table 5. Calculated values of X n for potential V 2 



n 


D 


d=2 


d=3 


d=5 


d=9 


2 


1 


.437 


.397 


.347 


.292 


4 


4 


.4428 


.4010 


.34912 


.29365 


6 


8,13 


.443007 


.401106 


.349171 


.293667 


8 


20,54 


.4430170 


.4011103 






10 


48 


.44301744 









We note how rapidly these numbers converge as one goes down each 
column in the tables. For each step increasing the order of approximation, 
we see one or two orders of magnitude increase in accuracy, somewhat better 
for E than for X. Also, one sees in these tables that the first approximation 
(a 'back of the envelope' computation) is accurate to about one percent. 
Such is the power of the X. For comparison, Table 6 presents results for the 
one-matrix problem, d=l and V\, computed by the same program. We see 
that the results of the many-matrix computations (above) converge about as 
rapidly as the one-matrix results, although the amount of work required to 
obtain the former is much greater. 



Table 6. Computed results for the one-matrix problem: V\ 



n 


D 


E 


X n 


2 


1 


.375 


.50 


4 


2 


.3717 


.5100 


6 


3 


.371638 


.51057 


8 


4 


.3716339 


.510611 


10 


5 


.37163373 


.5106136 
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Table 7 gives the E/d results computed for the potential V3 and one sees 
rapid convergence here as well. 



Table 7. Calculated values of E/d for potential V3 



n 


D 


d=2 


d=3 


d=5 


d=9 


2 


1 


.236 


.298 


.375 


.4725 


4 


4 


.2312 


.29470 


.373207 


.471358 


6 


8,13 


.231036 


.294625 


.3731823 


.47134965 


8 


20,54 


.2310258 


.29462242 






10 


48 


.23102504 









In Table 8 we see the results for the potential V4, which has the greatest 
amount of "flat directions" among our models. Here the rate of convergence 
is noticeably slower than in the previous models, but still looks convincingly 
good. 



Table 8. Calculated values of E/d for potential V4 



n 


D 


d=2 


d=3 


d=5 


d=9 


2 


1 


.24 


.30 


.38 


.47 


4 


4 


.224 


.289 


.370 


.4690 


6 


8,13 


.2232 


.2890 


.36944 


.468940 


8 


20,54,72 


.22299 


.28895 


.369431 




10 


48 


.222964 









Also, in the several tables above, one sees a suggestion of more rapid 
convergence for larger values of d; see the discussion of the large d expansion 
in Subsec. 5.3. 

In another experiment, we studied the one- matrix problem with potential 

V{<P) = \<P 2 - 9 -/ (7.1) 
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as the parameter g approached the value \/8/3ir where the bound state dis- 
appears. The numerical procedure searching to minimize the energy worked 
well until one approached very close to this critical value; then it failed dra- 
matically. 



Other X w values are also produced in these computations, albeit with a 
somewhat lesser accuracy. Table 9 has some of these for the potential V^. 



Table 9. Computed values of some other X w for V2 



X w 


d=2 


d=3 


d=5 


d=9 




-.0132659 


-.0082358 


-.004201 


-.001798 


X1122 


-.0066329 


-.0041179 


-.002101 


-.000899 


X-V2V1 


0.0 


0.0 


0.0 


0.0 



If there is rotational symmetry in the ground state, one can derive the fol- 
lowing relation among the fourth order X's, 

-Xllll — 2Xn22 + X\2V2 (7-2) 

and the data in Table 9 satisfy this relation, as does the corresponding data 
for the potential V4, which is also rotationally invariant. 



There is another, unexpected, phenomenon seen in the data of Table 9: 
namely that X1212 = 0. An increasing number of other X w 's also vanish when 
one looks at higher orders. This result also appears for the potential V3, but 
not for V4. When a particular X w goes to zero, so does the corresponding 
coefficient R w . The empirical rule is this: Write out the word w and remove 
any pair of matching adjacent letters; repeat this process; the X w will vanish 
unless this process can reduce the original word to null. I do not have a 
full explanation for this newly discovered selection rule but it appears to be 
related to the fact that these potentials (see (|6.2b|) and (|6.2c| )) involve only 
pairs (4>m4> m ) of each operator. This new symmetry is particular to Large N 
matrix mechanics with its noncommuting coordinate operators; it would not 
arise in ordinary quantum mechanics. 
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From an experimental (numerical) perspective, but lacking any formal 
proof, it appears that these types of large N problems are now solvable. It 
will be important for others to repeat this work independently in order to 
verify these results. 

8 Details of the Computational Program 
8.1 The Full F m 

The expression ( |6.3|) for F m ((f>) needs to be refined. The motivation for what 
follows comes from App. E in Ref. where the ground state wavefunction 
(the action) is modeled and one sees the consequent structure of F m (4>). 

Corresponding to each basic word Wi we want to have a group of terms 
(in the T w {4>)) with a common coefficient R^: 

F m (<P) = J2R ( r ) F m M- (8-1) 
For the first stage in this construction we define 

d m T w (<j))= J2 T -'(0) (8-2) 

w^mw' 

which, one can show, will guarantee that the flatness condition || 

[f m , F n {4>)\ - [vr n , Fjj>)\ = (8.3) 

is satisfied. 

For the second stage we take all permutations among the m = 1 . . . d 
letters that occur in the basic words Wj. 

F m .i{4>) = 7 w 1 , 7Tj d m Yl Permute T Wi (<P) (8.4) 
c[Wi) (a — 1)1 , 

v l > V / perm s 

where the constant c(w), the number of subcycles in the word w, is defined 
in App. A. The normalization constants used above are convenient but not 
essential. 

Now we construct the matrix elements 

(no sum) (8.5) 
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where these are linear combinations of the K WiV] i defined in ( |6.7| ) and Eq. ( |6.6| ) 
is replaced by 

Y^hjR^ = \ki> i = l...D. (8.6) 

3=1 

In order to save computing time in evaluating each tij it is important to find 
and to count repeated evaluations of the same K elements. I am not sure 
that I have done this job completely in my program. 



8.2 Searching 

The hardest part of this program is searching for the minimum energy in 
the parameter space of the basic word connected parts: X Wi = Xi, i = 
1 . . . D. The first method used fits a quadratic function to E(x) evaluated at 
D(D+l)/2 nearby points and then finds the extremum: 

bi = E(xi + 5) — E(x), Oij = E(x{ + 5, Xj + 8) — E(x) — hi — bj (8.7a) 

D 1 

a M v 3 =bi- -a^i, x[ = Xi- 8v { . (8.7b) 

3=1 1 

If one is close enough to the minimum, iterating this procedure should con- 
verge rapidly. For most of the data presented in Sec. 7 this method worked, 
although I am sure that more sophisticated techniques could have been more 
efficient. For the largest size computation carried out (d=5, n=8, D=72) 
the time for each evaluation of the energy was about one minute and each 
iteration of this search procedure took about 44 hours on a common desktop 
microcomputer. 

Sometimes, however, this approach failed. For the potential V4, beyond 
the sixth order calculation (for d=2 and d=3) this method diverged or led 
to impossible output (see the next subsection). What succeeded in those 
cases was a second method: start by solving the numerical problem for some 
other potential (like V3 where the first search method worked well) and then 
gradually change a coupling constant g inserted into the potential and solve 
again, repeating in small steps until one arrives at the desired result. At each 
new step one can start efficiently with a sort of perturbation theory 

f; a it jAj = 5 2 x[ = Xl - (5g)A i (8.8) 
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which involves the matrix a^j ( |8.7a ) which one has already calculated at the 
previous step. 

Just because the numerical search appears to converge is no proof that we 
have found the correct solution. In work on the potential V4 for d = 2 we had 
some results at the sixth order (D=8) which first appeared well converged 
by the first searching method; but a later check on the rotational symmetry 
( |7.2|) showed that this was a false solution. Repeating this calculation using 
the second search method described above led to satisfactory results. The 
fact that the false energy value was off only in the fifth decimal place stands 
as a cautionary note on this new numerical technique. 

Another numerical searching procedure is suggested by the algebraic work 
in Sec. 5. One could vary only the subset of X w 's with [w] = n*, keeping all 
others fixed; then cycle through the choices of n*. 

It should be repeated that this is all experimental work that is in need of 
sound mathematical justification and guidance. The multidimensional energy 
surface E(x) is a very complicated nonlinear function of the parameters x. 
In fact, there are singularities which may lie not far away from the desired 
minimum. One can see the simplest example of this situation in the 2x2 
matrix equation (|S.6|) for the d=l case. 



8.3 Constraints 

The quantities X w cannot be regarded as completely independent variables. 
For example, in the one-matrix case one has 

(0 I ((f) 2 - < <f) 2 >) 2 I 0) > (8.9) 

which leads to the inequality X 4 > — (X 2 ) 2 . 

Using the general Schwarz inequality, we can write 

|(0 I T W T W > I 0)| 2 < (0 I Tu,T w I 0)(0 I T-,T W , | 0) (8.10) 

for all words w and w' . This implies many constraints upon the allowed values 
of the X parameters as we search to minimize the energy. It is unclear how 
best to implement these constraints; in the computations reported here I 
only checked that the matrix ( |8.5| ) satisfied 

\ti,j\ 2 < u,i hv ki>° Vi,j (8.11) 
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at each evaluation. A failure of this test signals that the search has strayed 
into forbidden territory. 

An entirely different sort of constraint comes from the use of a purely 
real (rather than complex) representation for the operators. This implies 
that we should have X w = X* = X®. With the extensive symmetry of the 
problems studied here many of these constraints are automatic; but at the 
10th order for d = 2 and at the 8th order for d > 2, one finds some basic 
words that do not satisfy w ~ w. Rather than imposing this constraint, we 
are satisfied to find that this equality comes out in the numerical results. 

9 Excited States 

After the ground state problem is solved, we consider excited (adjoint) states 
in the reduced Hilbert space: 

H | 0) = E | 0), H\E) = E\ E) (9.1) 

where it should be remembered that we do not know the form of the reduced 
Hamiltonian H |J but only that it generates time translations. With the 
postulate 

I E) = U | 0) (9.2) 
for some operator U we find the identity 

(E - E Q )(0 I rfU | 0) = -i(0 | U ] U | 0). (9.3) 

Now we make the construction, as with F before, 

U = Y,rM<P), tf+ = ECW) ( 9 - 4 ) 

to w 

and we have, using ( |A.9|) , 

U = ^""^ t w T Wl TT m T W2 (9-5) 

w w=w\mv)2 

where the r w are as yet undetermined constants. 
We can now write ( |9.3| ) as 

E - E = ( Ab,«/?v)/ ( r w K w,w'r w ') (9.6) 

vi, w' w,w' 
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where the matrix K w ^ w i was defined earlier and from ( |A.10| ) we have 



L w>w i = — i(0 | T w t w i | 0) = - ^2 K U)V fK u r tV . (9.7) 

^ w=umv w '=u'mv' 



Finally, vary the coefficients r to find stationary values of ( |9.6| ) and we get 
a traditional linear matrix problem, where E — E is an eigenvalue of the 
matrix L with respect to the metric matrix K. 

The evaluation of the matrix K and thus also of L is done entirely in terms 
of the X w 's, which were already solved with the ground state problem. Thus 
(although I have not done any explicit numerical calculations for excited 
states) the complete spectrum of H can be calculated. The lowest order 
approximation, U = T m (0), gives E m - E = l/(2X mm ). 
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Appendix A. Useful Formulas Old and New 

Further conventions on the word notation: 
w = is the null word. 

w = m means that the word w consists of a single letter m. 
w ~ w' means that the two words differ by at most a cyclic permutation of 
their letters. 

w ~ w' means that the two words are equivalent under some larger symmetry. 
W1W2 = W3 means that the second word is appended to the first word and 

the result is the third word. 
w = umv means that the word w is decomposed as indicated. 
w is the word formed by reversing the sequence of letters in the word w. 
c(w), the number of subcycles in the word w, is defined as the largest integer 

k such that w = u k for any word u with [u] > 0. 

Basic relations among T(0) and X 0: 

T m w 0m-^ui ^ t X mWl T W2 (A. la) 

W=WlW2 
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T wm T w cj) m ^ ' T wl X W2 , m (A. lb) 

W = WlW2 

X mw = (0 I (j) m T w | 0) = (0 | T w (j) m | 0) = X wm (A.lc) 
Ti=T^ X* W = X W . (A.ld) 

Relation between X and Z: 

Z(j) = l + X(jZ(j)). (A.2) 

Examples (for the case of parity symmetry, which means that each letter 
must appear an even number of times or else the Z and X vanish): 

Zmn = (0 | m n I 0) = 5 mn X mm (A.3a) 
( X mmmm + 2X 2 mm if m = n = p = q 

^mnpq Zinpqrn \ X mm pp -\- X mm Xpp if 771 71 ^ p Q (A. 3b) 
[ Xmnrnn if p = 7JI ^ 71 = q. 

For one- matrix problems the label w is replaced by n = [w]. For systems 
with parity selection rule: 

T = l, T 1 = 0, T 2 = 2 -X 2 , X 2 =<0 2 > (A.4a) 
T 3 = 3 - 2<j)X 2 , T 4 = 4 - 30 2 X 2 - X A + X\ (A.4b) 
X 4 =< 4 > -2X 2 , X 6 =< 6 > -6X 4 X 2 - 5Xl (A.4c) 

Below are some new relations involving T(0) that are used in the present 
work. Start with the generating function 

Y = 1/(1 - [3 m( p m + X{(5)) = £ P w T w {<t>) (A.5) 

and calculate the commutator, 

i[n m ,Y] =Y$ m | 0)(0 | Y. (A.6) 
Now expand in powers of (3 and match terms to find: 

i[* m ,T w {4>)\ = ]T T W1 (0) | 0)(0 | 7U0). (A.7) 

ui=uiim«)2 

The other version of this relation 

i[n m ,f w ]= Yl f W2 \0)(0\f Wl (A.8) 

w=wimw2 



22 



comes from Eq.(D.ll) in Ref. 0. In a very similar way one gets the time 
derivative equation 

^-T w {4>) = T Wl n m T W2 (A. 9) 

where we have used ^4> m = 7T m . Combining the last two equations leads to 
i(0\T w ,jT w \Q) = ~ E (0 1^' | 0)(0 | T U ,T„ | 0) (A.10) 

w=umv w'=u'mv' 

which is surprisingly simple. 

Appendix B. Evaluating < T W T W > > 

We seek some recursive procedure for evaluation of the matrix elements 
K w>w , = (0 | T w {4>)T wl {<j)) | 0) = K w , >w (B.l) 



in terms of the connected parts X w . Using equations ( |A.la|) and ( |A.lb| ) it is 
relatively easy to find the following relations 

-K wm w ' K w mw ' -\- ^ X mu K w ^ v s X vrn K u ^ w i (B.2) 

w'=uv w=uv 

with the boundary counditions K Wt o = K 0tW = S Ws q. This looks very nice 
as a recursive computer program but it turns out to be expensive: the time 
required grows exponentially as one increases the size of the words. One 
could save time by building a table of all the K matrix elements one might 
need, but that requires enormous amounts of space. 

An alternative method is given by the following formula 

= E E Ku,v>X vu ,, [v] > 0, M > 0, K 0>0 = 1 (B.3) 



w=uv w'=u'v' 



which may be derived by combining equation ([A.la]) with the expansion 



YXrnviGM (B.4) 
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from Ref. J7| and also using the identity 

(0 | T W T W ,G^ | 0) = ]T (0 | T W1 T W2 | 0) 5 (B.5) 

w"=uv 

which is similar to the Ward Identities derived in App. E of Ref. 0. 

The program uses ( |B.3|) to build a small table of iT's each time one of 
them is called for and the time for this grows as n 4 rather than exponentially. 
Still, this is the main time consuming part of the computations. 



Appendix C. Some Alternative Computational Schemes 

One alternative scheme is to start out by fitting the quantity E mn (4>) 
instead of F m (4>), 

E mn (<P) = E Ri r ) T w (<P)- (CI) 

w 

The definition (|1.9d| ) is 

E mn {4>) I 0) = 2i[ff n , F m {(f))] | 0) (C.2) 
and using (|A.7| ), we find 

R ( r ] = 2R { 2 ] (G.3) 

upon comparison with (|6.3|) . Next, we use the formal expansion from Ref. [|7| 

(E l ) mn — y^,X wmn Gyj((f)) (C.4) 
w 

to write the system of conditions 

(0 | T w ,{E mn (E- l ) np - 5 m;P ) | 0) = (C.5a) 
X npw »(0 | T w 'T w Gu," | 0}R^ = -5m :P 5 w >fi (C.5b) 

w,w" 

and one can show, using (|B.3|), that this reduces to equations identical to 



(|6.6|) . So this method is not alternative at all. 

A second alternative scheme does away with minimizing the energy and 
works instead from the equations of motion: 

7r m | 0) = zF m (0) I 0) = -VM) I 0). (C.6) 
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Using the representation ( |6.3| ) for F m , this leads us to a new system of equa- 
tions 

- «5> I T w /f w | 0)^ = (0 | T W ,V^ | 0) (C.7) 

w 

where the matrix elements on the left side are the quantities L w i w defined 
in Sec. 9. One now has two sets of matrix equations - (|6.6|) and ( |C.7|) - 
determining the same set of expansion coefficients: call the solutions R and 
R' . One would now seek a set of values for the parameters X Wi that would 
make these two sets of solutions the same. Computationally, the way to do 
this would presumably be to minimize the error, 

Error =^| J R i -^| 2 (C.8) 

i 

and this defines another nonlinear search procedure. But what weight func- 
tion ought optimally to be put into this error calculation? 

A third alternative is to use the monomials <p w instead of the polynomials 
T w (4>) as a basis for the fitting of the operators F m or U. This leads to much 
simpler formulas for the matrix elements of K and L, expressed in terms of 
the moments Z w = (0 | (f) w | 0). Then one would use the relation (|A.2j) to 
evaluate each Z w in terms of the chosen set of parameters X w . I believe that 
this approach has drawbacks in both speed and numerical accuracy; but it 
should be explored. 



Appendix D. Is this Method Useful in Ordinary QM? 

With the apparent success of this approximation method in Large N 
matrix mechanics, one goes back to ordinary quantum mechanics to see if 
we have a new useful calculational technique. The formalism developed in 
Ref. is easily modified to fit the standard commutation relation 

i\Pi,Qj] = 8ijI (D- 1 ) 
with the following construction: 

y = e Pm-x(p) = J2 T^q) (D.2a) 

X{(3)=Y, C ^ X » ( D - 2b ) 
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Z(/3) = (0 I e Pm | 0) = J2 C ^(° I ^ I 0). (D.2c) 

Here \x represents the unordered set of occupation numbers {/x^} (remember 
that the g^'s commute with one another now) and 

C, = l/H(^)\. (D.3) 

i 

In the simple one-matrix case we have 

oo . 

= E P n Zn/n\, Z n = (0 | q n | 0) = / dq q n p(q) (D.4a) 

n=0 

oo 

X(P) = E P n XJn\ = ln(Z(P)) (D.4b) 

n=l 

and we want to test whether the ratio X n /Z n decreases rapidly with n, as we 
saw for the Large N situation in Sec. 2. For the case of a harmonic oscillator, 
we have the same result in both theories, namely X n vanishes for n > 2. 

One simple (non-oscillator) model that allows analytic calculations is a 
constant density p(q) over some finite range of q. Here we find that the ratio 
X n /Z n decays rapidly with n for the Large N situation but this ratio grows 
very rapidly for the ordinary quantum mechanics situation. 

We have also applied the method of this paper to the quantum mechanical 
nonlinear oscillator, 

H=\p 2 + \q A . (D.5) 

Numerical results for the ground state are shown in Table 10. The conver- 
gence seen here is fairly good, although not as good as for the similar Large 
N problem, shown in Table 6. (The accuracy shown here is comparable to 
that obtained with conventional variational calculations of this Schrodinger 
equation, at the same levels of approximation.) 



Table 10. Results for the Schrodinger equation (|D . 5| ) 



n 


D 


E 


X u 


2 


1 


.429 


.437 


4 


2 


.4217 


.4525 


6 


3 


.4210 


.45512 


8 


4 


.42086 


.45571 
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It must be reported, however, that the results shown in Table 10 were not 
obtained easily. The problem of nearby singularities in the energy surface, 
mentioned in Subsec. 8.2, was more severe in this ordinary quantum me- 
chanics problem than in the Large N problems. For the calculations through 
D = 3 I used the second searching method, starting from the harmonic os- 
cillator and then moving gradually to the quartic potential in steps of size 
1/8. For D = 4 I had to decrease the step size to 1/16; and for D = 5 I gave 
up after failing in the search procedure with step size 1/32. 

In conclusion, I am still in doubt about the answer to the question posed 
in the heading of this appendix. 
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